##########################################
########################################## 
## STUDY 1
##########################################
########################################## 
setwd(getwd())
source("libraries.R")

# load data for study
data_raw <- haven::read_dta(here("Data", "study1.dta"))

##########################################
########################################## 
## Recode Variables
##########################################
########################################## 
df <- data_raw %>% 
  rename(female = gender,
         age = age_0) %>% 
  mutate(race = case_when(
    race == 2 ~ 0,
    TRUE ~ 1
  ),
  
  #Party Identification (1 = Democrat; 0 = Republican)
  repdem = case_when(
    rep1 == 1 | rep1 == 2 | repdem2 == 1 ~ 0,
    dem1 == 1 | dem1 == 2 | repdem2 == 2 ~ 1,
    TRUE ~ NA_real_
  ),
  
  #Dummy indicating whether participant is Independent or not
  independents = case_when(
    repdem2 == 3 ~ 1,
    TRUE ~ NA_real_
  )) %>% 
  
  #Need for Chaos
  mutate(nfc = rowMeans(select(., c("needchaos_1", "needchaos_3", "needchaos_4", 
                                    "needchaos_5", "needchaos_6", "needchaos_9", "needchaos_10")), 
                        na.rm = TRUE)) %>% 
  
  #Political Violence Intentions
  mutate_at(vars(conflict2_2:conflict2_7), ~na_if(., 6)) %>% 
  mutate(violentact = rowMeans(select(., c("conflict2_2", "conflict2_3", "conflict2_5", "conflict2_6", "conflict2_7")), 
                               na.rm = TRUE)) %>%
  
  #Political Rumors: Believing and Sharing of DEMocratic and REPublican rumors
  mutate(true_DEM = rowMeans(select(., c("q15_1","q16_1","q17_1","q18_1","q19_1","q20_1")), na.rm = TRUE)) %>% 
  mutate(share_DEM = rowMeans(select(., c("q15_2","q16_2","q17_2","q18_2","q19_2","q20_2")), na.rm = TRUE)) %>% 
  mutate(true_REP = rowMeans(select(., c("q23_1","q24_1","q25_1","q26_1","q27_1","q28_1")), na.rm = TRUE)) %>% 
  mutate(share_REP = rowMeans(select(., c("q23_2","q24_2","q25_2","q26_2","q27_2","q28_2")), na.rm = TRUE)) %>% 
  
  #Standardize variables
  mutate(across(.cols = c(nfc, violentact, age, edu, true_DEM, share_DEM, true_REP, share_REP),
                ~ (scale(.) %>%  as.vector),
                .names = "{.col}_zscore")) %>% 
  
  #Scale variables to range from 0 to 1
  mutate(across(.cols = c(nfc, violentact, age, edu, true_DEM, share_DEM, true_REP, share_REP),
                ~ (zero1(.) %>% as.vector),
                .names = "{.col}_01"))

##########################################
########################################## 
## Sample Description 
## (See Supplementary Materials SM1, Table S2)
##########################################
########################################## 
#Gender
print(paste(round(mean(df$female), 2)*100-100 , "percent of participants are female"))
#Age
print(paste("Mean age of participants:" , round(mean(df$age) + 18, 2), 
            ", SD:", round(sd(df$age) , 2)))
#Education
print(paste("Less than high school:", round(prop.table(table(df$edu)), 2)[1]*100, "percent"))
print(paste("High school graduate:", round(prop.table(table(df$edu)), 2)[2]*100, "percent"))
print(paste(" Some college, but no degree:", round(prop.table(table(df$edu)), 2)[3]*100, "percent"))
print(paste("2 year college degree:", round(prop.table(table(df$edu)), 2)[4]*100, "percent"))
print(paste("4 year college degree:", round(prop.table(table(df$edu)), 2)[5]*100, "percent"   ))
print(paste("Graduate or professional degree:", round(prop.table(table(df$edu)), 2)[6]*100, "percent"))
#Ethnicity
print(paste("White/caucasian:", round(prop.table(table(df$race)), 2)[1]*100, "percent"))

##########################################
########################################## 
## Factor Analysis, etc.: Need for Chaos 
## (See Supplementary Materials SM2, Table S3)
##########################################
########################################## 
#Exploratory Factor Analysis: Need for Chaos - two-factor structure
df_EFA_two_factor <- df %>% 
  select(needchaos_1:needchaos_11) 

nfc.fa <- psych::principal(df_EFA_two_factor, 
                           nfactors = 2 , rotate="none" )  

#Display statistics for two factors
print("Factor loadings, two-factors:")
nfc.fa$loadings

#Exploratory Factor Analysis: Need for Chaos - final, one factor structure 
df_EFA_final_factor <- df %>% 
  select(needchaos_1, needchaos_3, needchaos_4, needchaos_5, needchaos_6, needchaos_9, needchaos_10) 

nfc.fa_final <- psych::principal(df_EFA_final_factor, 
                                 nfactors = 1 , rotate="none" )  

#Display statistics for final factor
print("Factor loadings, final factor:"); nfc.fa_final$loadings

#Cronbach's alpha for final factor
nfc <- df_EFA_final_factor %>% 
  psych::alpha(.) 
print(paste("Cronbach's alpha, Need for Chaos:" , round(nfc$total$raw_alpha,3)))  

##########################################
########################################## 
## Remove 5% with highest Need for Chaos?
## When set to TRUE, used for producing Fig S7, Fig S8, Fig S9
##########################################
##########################################
# Set to TRUE for excluding top 5% of participants with highest Need for Chaos
exclude_obs <- FALSE

if (exclude_obs == TRUE) {
  
  df <- df %>% 
    filter(nfc_zscore < quantile(nfc_zscore, probs = .95, na.rm = TRUE))
  
  print("Excluding top 5% observations with highest levels of the Need for Chaos")
  
}

##########################################
########################################## 
## Analysis - Fig. 2 of Main Text
##########################################
########################################## 
#Fig. 2: Left-Hand Panel
mean_nfc <- round(mean(df$nfc_01),2)
sd_nfc <- sd(df$nfc_01)

fig2_left <- ggplot(df, aes(x=nfc_01)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), color="grey", fill="grey") + 
  geom_vline(aes(xintercept=mean(nfc_01)),
             color="black",  size=1) +
  geom_vline(aes(xintercept=mean(nfc_01) + sd(nfc_01)),
             color="black",  size=1) +
  xlab("Need for Chaos") +
  ylab("Proportion") +
  theme(axis.text.x = element_text(size=13),
        axis.title = element_text(size=17),
        axis.text.y = element_text(size = 13)) +
  xlim(c(0,1)) + ylim(c(-.02,.15)) +
  annotate("text", x=mean_nfc + .075, y=-.015, label="Mean", color = "black", size = 5) +
  annotate("text", x= mean_nfc+sd_nfc + .075, y=-.015, label="+1 SD", color = "black", size = 5) 

#Save output
save(fig2_left,file=here("Data_recoded","fig2_left.R"))

##########################################
########################################## 
## Analysis - Fig. 3 of Main Text
##########################################
########################################## 
m_share_dem <- lm(share_DEM_01 ~ nfc_zscore + repdem  + female + edu_zscore + age_zscore + factor(race) , data = df )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "nfc_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study1")

m_share_rep <- lm(share_REP_01 ~ nfc_zscore + repdem  +  female + edu_zscore + age_zscore + factor(race) , data = df )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "nfc_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study1")

m_true_dem <- lm(true_DEM_01 ~ nfc_zscore + repdem +  female + edu_zscore + age_zscore + factor(race) , data = df )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "nfc_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study1")

m_true_rep <- lm(true_REP_01 ~ nfc_zscore + repdem + female + edu_zscore + age_zscore + factor(race) , data = df )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "nfc_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study1")

study1_coef <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study1_coef,file=here("Data_recoded","study1_fig3.R"))

#regression table for Supplemental Materials (Table S8)
stargazer(m_share_dem , m_share_rep , m_true_dem , m_true_rep , 
          align=TRUE,
          no.space=TRUE,
          column.labels=c("Share - Dem", "Share - Rep", "True - Dem", "True - Rep"),
          covariate.labels = c("Need for Chaos (std.)", "Party ID (1 = Dem.)", "Women", "Education", "Age  (std.)", 
                               "Ethnicity (1 = non-white)",  "Intercept"),
          keep.stat = c("n","ser"),
          dep.var.labels.include = FALSE)

##########################################
########################################## 
## Analysis - Fig. 4 of Main Text
##########################################
########################################## 
#Fig 4: Democratic Voters 
df_dem <- df %>% 
  filter(repdem == 1 )

m_share_dem <- lm(share_DEM_01 ~ nfc_zscore +  female + edu_zscore + age_zscore + factor(race) , data = df_dem )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study1")

m_share_rep <- lm(share_REP_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_dem )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study1")

m_true_dem <- lm(true_DEM_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_dem )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study1")

m_true_rep <- lm(true_REP_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_dem )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study1")

study1_coef_dem <- rbind(share_dem, share_rep, true_dem, true_rep)

#Fig 4: Republican Voters
df_rep <- df %>% 
  filter(repdem == 0 )

m_share_dem <- lm(share_DEM_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_rep )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study1")

m_share_rep <- lm(share_REP_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_rep )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study1")

m_true_dem <- lm(true_DEM_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_rep )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study1")

m_true_rep <- lm(true_REP_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_rep )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study1")

study1_coef_rep <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study1_coef_dem,file=here("Data_recoded","study1_fig4_dem.R"))
save(study1_coef_rep,file=here("Data_recoded","study1_fig4_rep.R"))

##########################################
########################################## 
## Analysis - Fig S2 of Supplementary Materials
##########################################
########################################## 
df_independents <- df %>% 
  filter(independents == 1 )

m_share_dem <- lm(share_DEM_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_independents )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study1")

m_share_rep <- lm(share_REP_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_independents )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study1")

m_true_dem <- lm(true_DEM_01 ~ nfc_zscore + female + edu_zscore + age_zscore + factor(race) , data = df_independents )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study1")

m_true_rep <- lm(true_REP_01 ~ nfc_zscore +  female + edu_zscore + age_zscore + factor(race) , data = df_independents )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "nfc_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study1")

study1_coef_independents <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study1_coef_independents,file=here("Data_recoded","study1_figS2.R"))

##########################################
########################################## 
## Analysis - Supplementary Materials SM8
##########################################
########################################## 
#Does Need for Chaos predict intentions to engage in political violence? 
m1 <- lm(violentact_01 ~ nfc_zscore + repdem  +  female + edu_zscore + age_zscore + factor(race), data = df)

#Table S27
stargazer(m1 , 
          align=TRUE,
          no.space=TRUE,
          column.labels=c("Political Violence Intentions"),
          covariate.labels = c("Need for Chaos (std.)" , "Party ID (1 = Dem.)" , "Women" , "Education" , "Age  (std.)" , "Ethnicity (1 = non-white)" ,  "Intercept"),
          keep.stat = c("n","ser"),
          dep.var.labels.include = FALSE)



